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We study numerically the phase-ordering kinetics of the site-diluted and bond-diluted Ising models 
after a quench from an infinite to a low temperature. We show that the speed of growth of the ordered 
domain’s size is non-monotonous with respect to the amount of dilution D: Starting from the pure 
case D = 0 the system slows down when dilution is added, as it is usually expected when disorder 
is introduced, but only up to a certain value D* beyond which the speed of growth raises again. 

We interpret this counterintuitive fact in a renormalization-group inspired framework, along the 
same lines proposed for the corresponding two-dimensional systems, where a similar pattern was 
observed. 


I. INTRODUCTION 

The quench of systems, like ferromagnets or binary mixtures, to below the critical point is characterized by the 
formation and subsequent growth of domains. The typical feature of the phase ordering process is dynamical scaling, 
whereby physical properties become time-independent if lengths are measured in units of L, i.e. the typical domain 
size at time t EHS- This scaling scenario is expected to be valid also in presence of quenched disorder when it does not 
prevent the phase-ordering process 0, [1|- In the case of pure systems power law growth L ^ is well established. 
Conversely, when disorder is present the nature of asymptotic growht has been much debated 

I n g eneral terms, disorder generated energy barriers bring about a slowing down of the coarsening process 
[T^ - fTsL [13, lUi - ll^l and, as a rule of thumb, more disorder means slower growth. However, recent studies have shown 
that the growth low does not behave monotonically with increasing disorder in the two dimensional site or bond 
diluted Ising model Q- In these systems, where a fraction D of sites or bonds is randomly removed from a regular 
lattice, one observes that, for sufficiently large D, more disorder produces faster growth. More precisely, starting from 
the pure case with D = 0, where the usual temperature-independent power-law L ~ is obeyed, one enters a region 
with an asymptotic logarithmic behavior as soon as D > 0. However a temperature-dependent power-law growth 
is observed right at D = Dc = 1 — Pc, where Pc is the critical percolation density, thus leading to a non-monotonic 
dependence of the speed of growth upon the disorder strength D. This counterintuitive behavior has been interpreted 
in terms of the topology of the network of occupied sites or bonds iO with the temperature-dependent algebraic 
behavior associated to the fractal properties of the percolating network. 

In this Article, by studying numerically the site-diluted Ising model (SDIM) and the bond diluted one (BDIM) in 
d = 3 we find a non-monotonous dependence of L on disorder qualitatively similar to the two-dimensional case. This 
suggests that the relation between growth-law and topology might be a rather generic property, a fact that could 
provide new insights into the problem of phase-separation in disordered ferromagnets. 

This paper is organized as follows. In Sec. |lT] we describe the models, the SDIM and the BDIM, that will be 
considered and studied throughout this paper. In Sec. IHII we present an overview of the results in d = 2 whereas 
numerical simulations of the d = 3 models are discussed in Sec. IlYl In Sec. El we conclude the paper with a final 
discussion of the results and of some open issues. 


II. THE DILUTED ISING MODEL 
A. The network 

In the SDIM a network is prepared by generating configurations of occupied sites as in random percolation. On 
the sites of a square three-dimensional lattice there are independent random variables n,, which take values = I 
(occupied site) with probability p and Ui = 0 (empty site) with probability D = 1 — p. The substrate is formed by 
the set of occupied sites. Any couple of nearest-neighbour occupied sites is connected by a bond. 

In the BDIM, instead, starting from a regular square three-dimensional lattice one removes bonds between nearest- 
neighbouring sites < ij > at random with probability D. Specifically, an independent random variable is defined 
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on each couple of nearest sites, which takes values Jij = 0 with probability D (bond is absent) and Jij = Jq with 
probability 1 — D (bond is present). 

For both models, D will be referred to as dilution. We now describe the geometrical structure of the network as 
D goes from high to low values. Recalling that the percolation threshold is Pc — 0.3116 and Pc — 0.2488 for site and 
bond percolation respectively, toi D > Dc = 1 — Pc, the network is a patchwork of finite clusters of occupied sites 
(bonds) with a characteristic size that diverges at the threshold and a fractal dimension df ~ 2.5 [l^. The dilution 
range D > Dc will not be considered in this paper since, lacking an infinite cluster, perpetual coarsening cannot be 
sustained. At Dc, the size of the percolating cluster of sites (bonds) diverges. As D is lowered below Dc, the infinite 
cluster has a fractal geometry over distances up to a typical length ^{D) that diverges at the threshold as 

i{D)^{Dc-D)-'', (1) 

with V ~ 0.875, whereas it becomes compact over larger distances. Moreover, clusters of finite size are also present. 
The structure of the network described above holds in the range Dc> D > D *, where D* is a dilution value where 

aD*)=a, (2) 

and a is a microscopic length, like the lattice spacing. For D < D*, the infinite cluster is compact on any lengthscale 
and finite clusters are absent. The remaining far-apart vacancies inside this cluster are essentially single-site objects. 
Their average distance defines a new characteristic length 

X{D) = a{D*/D)i. (3) 

B. Spin system 

The diluted Ising model is described by the following Hamiltonian 

H = - (4) 

<ij> 

Here Ui = ±1 are spins on the 3-d square lattice, and the couplings Kij for the SDIM is chosen as Kij = J^niUj, 
where Jq is a constant, and = 0 or = 1 with probability D and 1 — D, respectively. For the BDIM, instead, 
rij = 1 on all sites and Kij = 0 or Kij = Jo with probability D and 1 — D, respectively. Notice that, since Jq > 0, the 
model is ferromagnetic. 


C. Equilibrium states 

The equilibrium phase-diagram of the model is pictorially represented in Fig. [TJ For D < Dc, at low enough 
temperature the system exhibits ferromagnetic order. In the {D,T) plane there is a critical line Tc{D), which 
separates the paramagnetic from the ferromagnetic phase. The critical temperature Tc{D) , which in the following 
will be measured in units of ks/J, where ks is the Boltzmann constant, decreases from the pure value (Tc(0) ~ 4.5) 
as the dilution is increased and vanishes at Dc {Tc{Dc) = 0). 

D. Time Evolution 

Non-conserved dynamics[l|, 0 is implemented by evolving the spins with single-spin-flip transition rates of the 
Glauber form 

w{a, -cTj) = i [l - (Ti tanh(i7j^/T)] . (5) 

Here, is the local Weiss field obtained by the sum 

fff = E (6) 

KLi 


over the set of nearest-neighbors Li of i. 
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FIG. 1: (Color online) Pictorial representation of the equilibrium phase diagram of the diluted Ising model in d = 3. The bold 
blue line with circles and the bold indigo one with squares are the critical temperature Tc{D) of the SDIM and of the BDIM, 
respectively. The temperatures T = Tsimui where simulations will be performed are marked with horizontal dashed lines, blue 
with circles for the SDIM and indigo with squares for the BDIM. Vertical arrows represent the quenching processes considered 
in the simulations. 


The system is initially prepared in the infinite-temperature disordered state and, at the time t = 0, it is suddenly 
quenched to a finite temperature T. This temperature is set in the simulations to a value Tgi^nui = 1.1 for the SDIM 
and Tsimui = 0.75 for the BDIM. These values are chosen as a compromise: on the one hand we want to work in the 
limit of low-temperatures because the interplay between the two types of growth-law we are interested in is expected 
in the T ^ 0 limit Q. On the other hand very low temperatures are useless in simulations for at least a couple of 
reasons. The first is that at very low temperatures pinning effects make the dynamics very slow. As a result the 
growing length L{t,D) does not vary appreciably in the time range accessed by the simulations and it is difficult to 
discriminate between different growth-laws. The second related reason is that at low temperatures the pinning centres 
produce a stop and go behavior of interfaces which, in turn, gives rise to an oscillatory behavior on top of the growth 
of L(t,D). These oscillations prevent a reliable determination of the growth-law. A similar behavior is observed on 
fractal substrates [T^, [H, 

In Sec. m we will present the results of simulations of systems with different values of D in the range [0,IIc]. It 
must be noticed that, as discussed more detailedly in Q, for dilutions D = Dc or close to it the working temperature 
T = Tsimui is necessarily larger than Tc{D), see Fig. [TJ Hence we are quenching above the critical temperature 
Tc{Dc) = 0, the system will eventually relax to a disordered state with a finite spin coherence length ^a{T) (not to 
be confused with the substrate property ^{D)) and the growth of L{t,D) will eventually stop at t = teq{T) when 
L{t,D) ~ ^a{T). In a preasymptotic stage with L{t,D) ^a{T) coarsening takes place and, since ^a{T) diverges 
very fast for T —>■ 0, this can last for a very long time. The situation is similar to the one found in other ferromagnetic 
systems with Tc = 0 as, for instance, the one-dimensional Ising model [1^ We have checked that Tsimui is 
sufficiently low as to never observe the equilibration of the system in the range of simulated times, and that the 
system keeps on coarsening at any simulated time. 


E. Observables 

The observable quantity of interest in this paper is the typical domain size L{t,D). For a model defined on a 
disconnected network, as it happens for D > D*, phase ordering occurs independently on the various parts of the 
system and, correspondingly, different definitions of the growing length can be given. Following [l^. we determine 
this quantity from the inverse excess energy as 

L{t,D) = [E{t,D)-E^{D)]-\ (7) 

where E{t,D) = (H) is the energy at time t, the angular brackets denote a non-equilibrium ensemble average, taken 
over random initial conditions and over dynamical trajectories and disorder realizations, and Eao(D) is the energy 
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of the equilibrium state at the final temperature T. Eq. 0, which is often used to determine L{t,D) in non-diluted 
systems [l|, has, in the present case, the further advantage that the disconnected finite parts of the substrate which 
are already ordered do not contribute to the computation of L{t,D). Indeed each droplet is surrounded by empty 
sites in the SDIM or by null bonds in the BDIM and hence there is no excess energy associated with it when all spins 
are aligned, as discussed in I,®- Using Eq. 0, then, one computes L{t,D) only on the regions where growth is 
active, which is the kind of information we are interested in. 


III. GROWTH LAW 

A. Overview of the behaviour of the 2d-SDIM 

In I, [3 it was shown that the kinetic properties of the 2d-SDIM and BDIM are due to the presence of three fixed 
points, in the sense of the renormalization group, existing on the D axis. The first fixed point is the trivial one of the 
pure system, located at D = 0, and is associated to the usual power-law growth L{t, 0) ^ The second one is the 
percolative fixed point located at Dc and is characterized by a temperature-dependent power-law growth 

( 8 ) 

with C{T) > 2. This new power-law growth is due to the fact that diluted sites act as pinning centres where interfaces 
get stacked unless an activation energy AE is supplied by the thermal bath. According to the argument presented in 
E[i3 this energy scales as AE ^ AlnL(t, D). This is expected to be true for sufficiently low temperatures. Assuming 
an Arrhenius time t ~ exp[AE/(ksT)] to exceed such energetic barriers one obtains Eq. ([5]) with an exponent 

1/C(r) = kuT/A (9) 


proportional to the absolute temperature. 

The two fixed points discussed above are repulsive, in the sense that as soon as D ^ 0, or D ^ Dc the asymptotic 
dynamics is governed by a different fixed point, located at D = D* (roughly half way between D = 0 and D = Dc), to 
which a logarithmic increase of L{t, D*) is associated. Due to these fixed point structure, if the system is prepared with 
an intermediate dilution between D = 0 and D* (or equivalently between D* and Dc) a crossover pattern is observed 
at a certain time tcrossiT, D) from an initial transient regime governed by the nearest unstable fixed point (D = 0 or 
Dc respectively), with a power-law increase of the domains size, to a late regime controlled by the attractive point 
at D*, characterized by a logarithmic L{t,D). The crossover phenomenon occurs when L{t,D) reaches the typical 
size associated to dilution, which is either X{D) or ^{D) depending on 0 < D < D* or D* < D < Dc- Given that 
both A and ^ decrease approaching D* (Eqs. (11131) 1. the slowest possible growth, namely the one where tcross{T, D) is 
smaller, is obtained at D = D*. Therefore by comparing L(t, D) for different values of D one finds a non-monotonous 
behavior: The growth slows down in going from D = 0 to D = D* and then speeds up again when D is further 
increased from D = D* up to D = Dc- This can represent a practical tool to identify D*. 

In the next section we will show that a similar structure is observed also in Sd. 


IV. NUMERICAL RESULTS 

1. Simulation details 

The details of the simulations are as follows. We have considered a three-dimensional square lattice system of 
N = sites, with C = 240 for the SDIM and C = 150 for the BDIM. We have checked that, with this choice, no 
finite-size effects can be detected in the time-regime accessed by the simulations. For every choice of the parameters, 
we have performed a certain number (in the range 10-100) of independent runs with different initial conditions, 
disorder realizations and thermal histories in order to populate the non-equilibrium ensemble needed to extract 
average quantities. 

To speed up the simulations we have used a modified dynamics where spin flips in the bulk of domains, namely 
those aligned with all the nearest neighbors, are prevented. This modified dynamics does not alter the behavior of the 
quantities we are interested in, as has been tested in the 2d version of the models and in a large number of different 
cases We have checked that this is also true in the present study. 
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2. 3d results: SDIM 

Let us start by considering the SDIM. The time dependence of L{t,D), for various dilution values, is plotted in 
Fig. [2]for the 3d-SDIM. Notice that we plot L{t, D)/L{15, D) in order to normalize the curves at the microscopic time 
traicro — 15 when the late-stage scaling regime is entered. 

In the pure case the asymptotic power-law L{t, 0) ^ is expected. In the relatively small range of simulated 
times we measure a somewhat smaller exponent (in the last decade a best fit procedure yields I/z ~ 0.46, as already 
reported in [^, , because the convergence to the true asymptotic behavior is very slow, as explained in [s^ . Indeed 

the observed effective exponent 1/z appears to increase with time. However here we are mainly interested in the 
non-monotonic behavior of the speed of growth for which an accurate estimate of z is not necessary. In fact, it is 
clearly observed that, by increasing D, the growth law becomes slower and slower up to a certain value of D cz 0.5 
that, following Q, we identify with D*. For D > D* growth becomes faster upon increasing D up to Dc- This 
behavior results in a non-monotonic D-dependence of L(t,D), at fixed time t. This can be clearly appreciated in 
the inset of Fig. where the behavior of the effective exponent obtained by fitting the data for various D with a 
power-law (regardless of the fact that the actual growth-law might not be power-law) in the time-range t > 5000. 

Having established the non-monotonous dependence of the speed of growth on D, we make some comments on the 
form of the growth-law. At the dilution D = D* we observe a downwards banding of L{t, D) which indicates that the 
growth is slower than algebraic, suggesting in this case a logarithmic law. We recall that a form L{t,D) ~ 
was predicted in a model with bond disorder in plj |. which, however, due its slow character has never been accurately 
demonstrated, to the best of our knowledge, even in the less numerically demanding 2d cases. Clearly in this 3d 
case the task is much more demanding, and the time range accessible in the simulations is not sufficient for a precise 
determination of the analytical form of the growth-law. 

For D = Dc, on the other hand, the behavior of L{t, Dc) is well fitted, as expected, by a power-law L{t, Dc) — ^ 

with l/C,{Tsimui) — 0.41. According to the discussion of Sec. IHI A[ this behavior is expected irrespective of the 
quenching temperature (but with a different temperature-dependent exponent Cfj, provided it is sufficiently low. In 
order to check this we have performed a series of simulations by fixing D = Dc and changing Tgimui- The behavior 
of L{t, D) in these cases is shown in Fig. [31 The data are consistent with power-laws for any choice of Tgimui- Notice 
however that for very small Tsimui the power-law is decorated by a sort of log-time periodic oscillation. The exponent 
1/C(T), obtained by fitting all the curves in the last decade (for T = 0.5 we use two decades in order to smear out the 
oscillation), is plotted in the inset. According to Eq. ([9|), one finds a linear behavior at sufficiently low temperatures. 
Fitting this curve we obtain A ~ 0.36 (Eq.®). Notice that for large values of the temperature this linear behavior 
breaks down, due to the constraint l/(f{T) <1/2 since it is impossible to grow faster than in the pure case. 


3. 3d results: BDIM 

In order to check the generality of the behaviors observed in the SDIM we have performed a series of simulations 
also for the 3d-BDIM. Unfortunately, for this model it is not possible to find a region of sufficiently low temperatures 
to be considered representative of the T —0 limit, as discussed in Sec. IHDl but high enough to eliminate the stop 
and go behavior of the interfaces, which is responsible for the oscillatory character of L(t,D). We have explored 
several values of Tgimui < 1, and in all cases we observe a non-monotonic behavior of L(t,D) as D is varied. In Fig. 
dH) we plot results for Tgimui = 0.75. For quenches at this temperature oscillations are clearly evident ruling out the 
possibility to make any statement on the form of the growth-law (whether logarithmic or power-law). Despite this, the 
non monotonic behavior of L{t, D) as D is varied at fixed t is quite evident from Fig. jj] Starting from the pure case 
D = 0, where a fast growth is observed (although much slower than the asymptotic expected behavior L(t,0) ^ 
due to the long-lasting pre-asymptotic corrections discussed above) the growth gets slower as D is increased up to 
D = D* ~ 0.5, and then speeds up again until Dc is reached. This shows that the non-monotonicity of the speed of 
growth as disorder is increased is a general property of diluted system, irrespective of dimension and of the kind of 
dilution. This suggests that an interpretation based on the presence of three fixed points, the pure one and those at 
D* and D^ and the attractive/repulsive nature of such fixed points, can be appropriated quite generally. 


V. CONCLUSIONS 

In this paper we have studied the diluted Ising model in three dimensions. The aim of the paper was to show the 
gnerality of the phase-ordering behavior previously observed in the 2d-SDIM and 2d-BDIM, where it was first shown 
I [a the non-monotonous dependence of the speed of growth on the amount of disorder, a fact that was interpreted as 
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FIG. 2: (Color online) The behavior of the typical length L{t, D)/L{15, D) after a quench to T — Tsimui is shown for the SDIM 
for several dilutions D (see legend) in a double logarithmic plot. The dashed-black line is the pure growth-law L(t,0) ^ 
and the dashed-violet one is the power-law L{t,Dc) ^ In the inset the effective exponent 1/zeff (see text) is plotted 

against D. 


due to the close relation between the coarsening properties and some topological features of the network in disordered 
magnets. In order to do that we have considered three-dimensional systems. 

Our results show that the pattern of behavior observed in these models is analogous to the one found in the 2d 
SDIM and RBIM, and henceforth an analogous interpretation can be given. In the pure system phase-ordering occurs 
without activation energy by a smooth curvature-driven displacement of interfaces, leading to the well known growth- 
law L{t,0) ^ The situation changes radically as soon as disorder is introduced because inhomogeneities pin 

the interfaces in energy minima and thermal activation is required to proceed. In [9| an argument was developed 
leading to the conjecture that the L{t, D)-dependence of the pinning barriers changes when passing from D < Dc 
to D = Dc, due to a radically different large-scale topology of the network which turns from a compact to a fractal 
object. This implies also a different growth-law, logarithmic for D < Dc vs power-law L{t, Dc) ~ right at Dc- 

In a renormalization-group language this features can be accomodated in a framework with three fixed points, located 
a.t D = 0, D = D* and D = Dc- Accordingly, a crossover phenomenon is observed where, for any 0 < D < Dc, 
L(t, D) changes from an early behavior associated to the nearest repulsive fixed point to the asymptotic logarithmic 
law. This gives rise to the non-monotonic behavior of L{t,D) as D is varied keeping t fixed. 

The presented numerical study suggests that the pattern of behaviors observed in is quite a general featnre 

of ferromagnetic disordered systems. Indeed we find that the non-monotonic behavior of L(t,D) is observed also 
in Sd systems. Our results also indicate that, at least for the SDIM, L{t, Dc) exhibits a power law growth with a 
growth exponent (^{t) inversely proportional to the temperature. The growth law becomes slower as D is decreased 
to D*, even if our data do not allows the precise determination of a logarithmic growth for L{t, D*). It would be an 
interesting line of research to investigate if a similar structnre underlies also the kinetics of other disordered systems 
snch as, for instance, Ising models with random fields. 
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